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1 DG discretization of the linear model problem 


Many engineering problems such as chemical reaction processes, heat conduction, 
nuclear reactors, population dynamics etc. are governed by convection-diffusion- 
reaction partial differential equations (PDEs). The general model problem used in 
the code is 



(la) 

(lb) 

(l c) 


The domain Q. is bounded, open, convex in with boundary dO. = U F^, 
F^ n F^ = 0, 0 < £ <C 1 is the diffusivity constant, f G L?' (fl) is the source func¬ 
tion, b G (lF^'”(fl))^ is the velocity field, is the Dirichlet boundary 

condition, ^ G //^/^(F^) is the Neumann boundary condition and n denote the 
unit outward normal vector to the boundary. 

The weak formulation of ([Til reads as: find uGU such that 



where the solution space U and the test function space V are given by 

U = {uGH\D.) : u = g^onr^}, V = {v G H\D.) : v = 0onF°}. 

Being the next step, an approximation to the problem ^ is found in a finite¬ 
dimensional space V/,. In case of classical (continuous) FEMs, the space Vh is set 
to be the set of piecewise continuous polynomials vanishing on the boundary BQ,. 
In contrast to the continuous FEMs, the DGFEMs uses the set of piecewise poly¬ 
nomials that are fully discontinuous at the interfaces. In this way, the DGFEMs 
approximation allows to capture the sharp gradients or singularities that affect the 
numerical solution locally. Moreover, the functions in F/, do not need to vanish at 
the boundary since the boundary conditions in DGFEMs are imposed weakly. 

In our code, the discretization of the problem ([Til is based on the discontinuous 
Galerkin methods for the diffusion part HJlhl and the upwinding for the convection 
part Em. Eet {(§/,} be a family of shape regular meshes with the elements (tri¬ 
angles) Ki G satisfying Q. = UA' and KinKj = 0 for K,, Kj G ^h- Let us denote 
by Fq, F/) and F^^ the set of interior, Dirichlet boundary and Neumann boundary 
edges, respectively, so that FoUF/jUF^? forms the skeleton of the mesh. For any 
K G^h, let ^k{K) be the set of all polynomials of degree at most k on K. Then, set 
the finite dimensional solution and test function space by 


Vh = {vG L\a) : v|j, G P^(^), G V. 
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Figure 1: Two adjacent elements sharing an edge (left); an element near to domain 
boundary (right) 


Note that the trial and test function spaces are the same because the boundary con¬ 
ditions in discontinuous Galerkin methods are imposed in a weak manner. Since 
the functions in Vh may have discontinuities along the inter-element boundaries, 
along an interior edge, there would be two different traces from the adjacent ele¬ 
ments sharing that edge. In the light of this fact, let us first introduce some notations 
before giving the DG formulation. Let Ki, Kj G i§/, (/ < j) be two adjacent elements 
sharing an interior edge e = KiCiKj cTo (see Fig[T]). Denote the trace of a scalar 
function v from inside K, by v, and from inside Kj by vj. Then, set the jump and 
average values of v on the edge e 

[v] = Ville - V jlle, {v} = ^ {Vi + Vj), 

where is the unit normal to the edge e oriented from Ki to Kj. Similarly, we set 
the jump and average values of a vector valued function q on e 

[q] = q,- • - qy • n^, {q} = ^(q; + qy), 

Observe that [v] is a vector for a scalar function v, while, [q] is scalar for a vector 
valued function q. On the other hands, along any boundary edge e = KjCi dO., we 
set 

[v] = v,-n, {v} = Vi, [q] = q,- • n, {q} = q,- 

where n is the unit outward normal to the boundary at e. 

We also introduce the inflow parts of the domain boundary and the boundary of a 
mesh element K, respectively 

F^ = {.r G dD. : b(.x) • n(.r) < 0} , dK = {x & dK : b(.r) • nK{x) < 0}. 

Then, the DG discretized system to the problem ([Til combining with the upwind 
discretization for the convection part reads as: find Uh G V/, such that 

ah{uh,Vh) = lh{vh) Vv/, G Vh, (3) 
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‘^h{uh,Vh) = / £'Vuh-'Vvhdx+ 52 / ■yuh +0CUh)vhdx 

Ke^i, Jk Ke^,^ Jk 

- 52 {£'^Uh}-[vh]ds+K Y, {syvh}-[uh]ds 

e^VoUT'o’'^ e^VoUl'i)^ 

+ 52 /" h-n{ul''‘— u'l^)vhds— 52 [ h-nu^j^Vhds 

+ 52 ^ /[m/i] • 

eerouro 

lh{vh)= 52 / fvhdx+ 52 i^Vh-eVvh-n]ds 
K%Jk e&oJe ) 

- 52 / b • + 52 / g^nds, 

K€^,,ddK-r\r- eeT^Je 

where and denotes the values on an edge from outside and inside of an el¬ 
ement K, respectively. The parameter fc determines the type of DG method, which 
takes the values {—1,1,0}: K = —I gives "symmetric interior penalty Galerkin" 
(SIPG) method, fC = 1 gives "non-symmetric interior penalty Galerkin" (NIPG) 
method and tc = 0 gives "inconsistent interior penalty Galerkin " (IIPG) method. 
The parameter a G Mq is called the penalty parameter which should be sufficiently 
large; independent of the mesh size h and the diffusion coefficient e [6i [Sec. 
2.7.1]. In our code, we choose the penalty parameter a on interior edges depend¬ 
ing on the polynomial degree k as a = 3k{k + 1) for the SIPG and IIPG methods, 
whereas, we take a = 1 for the NIPG method. On boundary edges, we take the 
penalty parameter as twice of the penalty parameter on interior edges. 

2 Descriptions of the MATLAB code 

The given codes are mostly self-explanatory with comments to explain what each 
section of the code does. In this section, we give the description of our main code. 
The use of the main code consists of three parts 

1. Mesh generation, 

2. Entry of user defined quantifies (boundary conditions, order of basis efc.), 

3. Forming and solving fhe linear sysfems, 

4. Plotting fhe solutions. 

Excepf fhe lasf one, all fhe parfs above fake place in fhe m-file Main_Linear.m 
which is fhe main code fo be used by fhe users for linear problems wifhouf need 
fo enfry fo any ofher m-file. The lasf pari, plotting fhe solufions, lakes place in fhe 
m-file dg_error.m. 
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2.1 Mesh generation 


In this section, we define the data structure of a triangular mesh on a polygonal 
domain in The data structure presented here is based on simple arrays [T] 
which are stored in a MATLAB "struct" that collects two or more data fields in 
one objecf fhaf can fhen be passed fo roufines. To obfain an inifial mesh, firsfly, 
we define fhe nodes, elemenfs, Dirichlef and Neumann conditions in fhe m-file 
Main_Linear.m, and we call fhe getmesh function fo form fhe inifial mesh sfrucfure 
mesh. 

% Generate the mesh 
% Nodes 

Nodes = [0 ,0;0.5,0;1 ,0;0,0.5;0.5 ,0.5;1 ,0.5;0,1;0.5 ,1;1 ,1]; 

% Elements 

Elemenfs = [4,1 ,5 ; 1 ,2 ,5 ;5 ,2,6; 2,3 ,6 ;7 ,4 ,8 ;4,5 ,8 ; 8 ,5 ,9;5 ,6 ,9]; 
% Dirichlet bdry edges 

Dirichlef = [ 1 ,2 ;2,3 ; 1 ,4; 3 ,6 ;4 ,7 ;6 ,9 ;7 ,8 ; 8 ,9]; 

% Neumann bdry edges 
Neumann = [ ] ; 

% Initial mesh struct 

mesh = gefmesh (Nodes , Elemenfs , Dirichlef ,Neumann); 

As if can be understood, each row in fhe Nodes array corresponds fo a mesh node 
wifh fhe firsf column keeps fhe x-coordinafe of fhe node and fhe second is for fhe 
y-coordinafe, and fhe i — th row of fhe Nodes array is called fhe node having index 
i. In fhe Elements array, each row wifh 3 columns corresponds fo a friangular ele- 
menf in fhe mesh confaining fhe indices of fhe nodes forming fhe 3 verfices of fhe 
friangles in fhe counfer-clockwise orienfafion. Einally, in fhe Dirichlet and Neu¬ 
mann arrays, each row wifh 2 columns corresponds fo a Dirichlef and Neumann 
boundary edge confaining fhe indices of fhe sfarfing and ending nodes, respectively 
(see Eigllll. The mesh "sfrucf" in fhe code has fhe following fields: 

• Nodes, Elemenfs, Edges, infEdges, DbdEdges, NbdEdges, infEdges 

• verficesl, verfices2, verfices3, 

• Dirichlef, Neumann, EdgeEls, ElemenfsE. 

which can be reached by mesh.Nodes, mesh.Elements and so on, and fhey are used 
by fhe ofher funcfions fo form fhe DG consfrucfion. The inifial mesh can be uni¬ 
formly refined several fimes in a "for loop" by calling fhe function uniformrefine. 

for jj=l:2 

mesh= uniformrefine (mesh ); %Refine mesh 

end 
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Figure 2: Initial mesh on the unit square Q. = [0,1]^ with nodes n,-, triangles Ej and 
edges ek 

2.2 User defined quantities 

There are certain input values that have to be supplied by the user. Here, we will 
describe that how one can define these quantities in the main code Main_Linear.m. 

One determines the type of the DG method (SIPG, NIPG or IIPG) and the order of 
the polynomial basis to be used by the variables method and degree, respectively. 
According to these choices, the values of the penalty parameter and the parame¬ 
ter {—1,1,0} defining DG method in ([3ll are set by calling the sub-function 
set_parameter. 

% method : NIPG = 1, SIPG = 2, 11PG=3 
method =2; 

% Degree of polynomials 
degree = 1; 

% Set up the problem 

[penalty , kappa]= set_parameter (method , degree ); 

The next step is to supply the problem parameters. The diffusion constant £, the 
advection vector b and the linear reaction term a are defined via the sub-functions 
fdiff,fadv and/react, respectively. 
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%% Define diffusion , advection , and reaction as subfunctions 


% D iffu si on 

function diff = fdiff(x,y) 

diff = (10^( — 6)). * ones ( size (x) ) ; 
end 

% Adv ection 

function [advl,adv2] = fadv(x,y) 
advl =( 1/ sqrt(5)) * ones( size (x)); 
adv2 =(2/ sqrt(5)) * ones( size (x ) ); 
end 

% Linear reaction 
function react = freact(x,y) 
react = ones( size( x )); 
end 

The exact solution (if exists) and the source function / are defined via the sub¬ 
functions fexact and fsource, respectively. Finally, the boundary conditions are 
supplied via the sub-functions DBCexact and NBCexact. 


% First derivative wrt x 

yex_x = ( — l./(sqrt(5*diff))).*(sech ((2* x—y —0.25)./... 

(sqrt(5*diff)))).^2; 

% First derivative wrt y 

yex_y = ((0.5)./(sqrt(5*diff))).*(sech ((2* x—y —0.25)./... 

(sqrt(5*diff)))).^2; 

% Second derivative wrt x 

yex_xx = ((0.8)./diff).* tanh ((2* x—y —0.25)./(sqrt(5*diff))).*... 

(sech((2* X—y — 0.25)./(sqrt(5*diff)))).^2; 

% Second derivative wrt y 

yex_yy = ((0.2)./diff).* tanh ((2* x—y —0.25)./(sqrt(5*diff ))).*... 
( sech ((2 * X—y — 0.25)./(sqrt(5*diff)))).^2; 

% Force function 

source=— diff . * ( yex_xx-i-yex_yy )-i-( adv 1 . * yex_x-i-adv2 . * yex_y ) -i-. .. 

2.3 Forming and solving linear systems 

To form the linear systems, firstly, let us rewrite the discrete DG scheme (l3]l as 
ah{uh,Vh) := Dh{uh,Vh) + Ch{uh,Vh) + Rh{uh,Vh) = h{vh) Vv/, € F/,, (5) 
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where the forms Dh{uh,Vh), Ch{uh,Vh) and Rf,{uh,Vh) corresponding to the diffu¬ 
sion, convection and linear reaction parts of the problem, respectively 


Dh{uh,Vh)=Y, £yuh-Vvhdx+ Y, ^ [uh]-[vh]ds 

- Y -[vhjds+K Y {£yvh}-[uh]ds 

Ch{uh,Vh)=Y ^■'^UhVhdx 

+ Y [ Y [ b- 

Ke^lrddK \da K€^,,ddK flF 

Rh{uh,Vh) =Y Cf'UhVhdx 

h{vh) =Y f^hdx -h y ( ^Vh - eVvh ■ n ) di 

K%Jk e&oJe \he ) 

Y [ ■ ^g^nds + Y [g^^hds, 

K€^i, d dK- nr- epr« de 


nu'^Vhds 


For a set of basis functions spanning the space Vh, the discrete solution 

Uh € Vh is of the form 

N 

Uh = Yj "^j^j ^d) 

i=i 

where u = (ui,U 2 , • • •, is the unknown coefficient vector. After substituting 
(iTl) into (l5]l and taking Vh = <pi, we get the linear system of equations 

N N N 

Y^ 0,) + Y^ 0!') + Y^ ; i — 1)2,.. . ,N 

j=i ;=i j=i 

( 8 ) 

Thus, for / = 1,2,... ,A, to form the linear system in matrix-vector form, we need 
the matrices D,C,R € related to the terms including the forms Dh, Ch and Rh 
in ([8]l, respectively, satisfying 


Dv+Cv +Rv = F 

with the unknown coefficient vector u and the vector F G related to the linear 
rhs functionals lh{^i) such thati^- = Ih{<pi), / = 1,2,... ,N. In the code Main_Linearm, 
all the matrices D,C,R and the vector F are obtained by calling the function global_system, 
in which the sub-functions introduced in the previous subsection are used. We set 
the stiffness matrix. Stiff, as the sum of the obtained matrices and we solve the 
linear system for the unknown coefficient vector coef:= u. 


%Compute global matrices and rhs global vector 
[D,C,R,F] = global_sy stem (mesh , @fdiff , @fadv , @freact , . . . 

@fsource , @DBCexact, @NBCexact, penalty , kappa , degree ) 

Stiff=D+C+R; % Stiffness matrix 

coef = S tif f \F; % Solve the linear system 

2.4 Plotting the solution 

After solving the problem for the unknown coefficient vector, the solutions are 
plotted via the the function dg_error, and also the L^-error between the exact and 
numerical solution is computed. 

% Compute L2—error and plot the solution 
[ 12err ,hmax]= dg_error(coef ,mesh , @fexact , @fdiff , degree ); 


3 Models with non-linear reaction mechanisms 

Most of the problems include non-linear source or sink terms. The general model 
problem in this case is 


au — eAu -|- b • Vm -|- r{u) 

= / 

in fl. 

(9a) 

u 


on F^, 

(9b) 

eVu ■ n 

= / 

on F^. 

(9c) 


which arises from the time discretization of the time-dependent non-linear diffusion- 
convection-reaction equations. Here, the coefficient of the linear reaction term, 
a > 0, stand for the temporal discretization, corresponding to 1 /At, where At is the 
discrete time-step. The model (|9ll differs from the model ([B by the additional non¬ 
linear term r{u). To have a unique solution, in addition to the assumptions given 
in Section 1, we assume that the non-linear reaction term, r{u), is bounded, locally 
Lipschitz continuous and monotone, i.e. satisfies for any > 0, 5 , 51,52 € K 

fhe following conditions [71 


|r;(5 )|<C, 00 

Ik,-(5i)-ri(52)||i2(ti) <F||5 i -52||i2(j2), L>0 

r,-€C'(]R+), r,-(0)=0, r;(5)>0. 

The non-linear reacfion ferm r{u) occur in chemical engineering usually in the form 
of products and rational functions of concentrations, or exponential functions of the 
temperature, expressed by the Arrhenius law. Such models describe chemical pro¬ 
cesses and they are strongly coupled as an inaccuracy in one unknown affects all 
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the others. 


To solve the non-linear problems, we use the m-file Main_Nonlinear which is 
similar to the m-file Main_Linear, but now we use Newton iteration to solve for 
/ = 1,2 ,... the non-linear system of equations 

N N N r 

ttyZ)/;( VjCh(^(j)j, (jii')Vj^ (^i)-\- I r(^Ufi)(j)idx — (10) 

j=l ' j=l j=l Jo. 

Similar to the linear case, the above system leads to the matrix-vector form 

Du-l-Cu-|-/?u-|-//(u) =F 

where, in addition to the matrices D,C,R ^ and the vector F G R^, we also 

need the vector H G R^ related to the non-linear term such that 

Hi{v) = / M L Mj 
Jo \j=l 

We solve the nonlinear system by Newton method. For an initial guess u° = 
(u°, U 2 , •.., we solve the system 

( 11 ) 

= >v*^ + u^ A: = 0,l,2,... 

until a user defined folerance is safisfied. In (fTTl) . Res^ and J^ denofe fhe vector of 
sysfem residual and ifs Jacobian mafrix af fhe currenf iferafe u^, respecfively, given 
by 


^ (j)idx , / = 1,2,... ,N. 


Res’^ = {D + C + R)v’^ + H{v’^)-F 
J^ = D + C + R + HJ{v'^) 


where HJ{v^) is fhe Jacobian mafrix of fhe non-linear vector H af 


HJ{v^) = 


5U|(te) 




(5i)| 






(5i)| 

. 


In fhe code Main_Nonlinear, obfaining fhe mafrices D,C,R and fhe rhs vecfor F is 
similar fo fhe linear case, buf now, additionally, we infroduce an initial guess for 
Newton iferafion, and we solve fhe nonlinear system by Newfon mefhod. 
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% Initial guess for Newton iteration 
coef=zeros(size(Stiff ,1) ,1); 

% Newton iteration 

noi=0; 

for ii =1:50 

noi = noi +1; 

% Compute the nonlinear vector and its Jacobian matrix at 
% the current iterate 

[H,HJ]=nonlinear_global (coef , mesh ,@freact_nonlinear , degree); 

% Form the residual of the system 
Res = Stiff*eoef + H — F; 

% Form the Jacobian matrix of the system 
% (w.r.t. unknown coefficients coef) 

J = Stiff + HJ ; 

% Solve the linear system for the correction "w" 
w = J \ (—Res); 

% Update the iterate 
eoef = eoef + w; 

% Check the accuracy 
if norm( J*w+Res) < le—20 

break; 
end 

end 

To obtain the non-linear veetor F[ and its Jaeobian F[J at the eurrent iterate, we 
call the function nonlinear_global, and it uses the function handle freact_nonlinear 
which is a sub-function in the file Main_Nonlinear. The sub-function/reacf_non/mear 
has to be supplied by user as the non-linear term r{u) and its derivative /(«). 

% Nonlinear reaction 

function [r,dr] = freact_nonlinear (u) 

% Value of the nonlinear reaction term at the current iterate 
r = u.^2; 

% Value of the derivative of the nonlinear reaction 
% term at the current iterate 
dr = 2*u; 

end 
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4 MATLAB routines for main code 


Here, we give the main m-file Main_Nonlinearm of the eode. The full eode is 
available upon request to the e-mail address uzuncaOgmail. com. 

1 % This routine solves the diffusion —convection —reaction equation 

2 % 

3 % \alpha u — \ epsilon *\Delta u + b\dot\nabla u + r(u) = f 

4 % 

5 % using DG—FEM. 

6 

7 function Main_Nonlinear () 

8 

9 clear all 

10 clc 

11 

12 % Generate the mesh 

13 

14 % Nodes 

15 Nodes = [0 ,0;0.5,0;1 ,0;0,0.5;0.5 ,0.5;1 ,0.5;0,1;0.5 ,1; 1 ,1]; 

16 % Elements 

17 Elements = [4,1 ,5 ; 1 ,2 ,5 ;5 ,2 ,6; 2 ,3 ,6;7 ,4,8 ;4 ,5 ,8 ; 8 ,5 ,9;5 ,6 ,9]; 

18 % Dirichlet bdry edges 

19 Diriehlet = [ 1 ,2;2 ,3 ; 1 ,4 ;3 ,6;4 ,7 ;6 ,9 ;7 ,8 ; 8 ,9]; 

20 % Neumann bdry edges 

21 Neumann = []; 

22 % Initial mesh struct 

23 mesh = getmesh (Nodes , Elements , Diriehlet ,Neumann); 

24 

25 for jj =1:2 

26 mesh=uniformrefine (mesh ); %Refine mesh 

27 end 

28 

29 % method : NIPG = 1, SIPG=2, IIPG=3 

30 method = 2; 

31 

32 % Degree of polynomials 

33 degree = l; 

34 

35 % Set up the problem 

36 [penalty , kappa]= set_parameter (method , degree ); 

37 

38 %Compute global matrices and rhs global vector 

39 [D,C,R,E]= global_system (mesh , @fdiff , @fadv , @freaet ,... 
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40 @fsource , @DBCexact, @NBCexact, penalty , kappa , degree ); 

41 

42 Stiff=D+C+R; % Stiffness matrix 

43 

44 % Initial guess for Newton iteration 

45 eoef=zeros ( size ( S tiff , 1), 1); 

46 

47 % Newton iteration 

48 noi=0; 

49 for ii =1:50 

50 noi = noi + l; 

51 

52 % Compute the nonlinear vector and its Jacobian matrix at 

53 % the current iterate 

54 [H,HJ]=nonlinear_global(eoef, mesh ,@freaet_nonlinear , degree); 

55 

56 % Form the residual of the system 

57 Res = Stiff*eoef + H — F; 

58 

59 % Form the Jacobian matrix of the system 

60 % (w.r.t. unknown coefficients coef) 

61 J = Stiff + HJ ; 

62 

63 % Solve the linear system for the correction "w” 

64 w = J \ (—Res ); 

65 

66 % Update the iterate 

61 eoef = eoef + w; 

68 

69 % Check the accuracy 

70 if norm ( J=i<w+Res ) < le—20 

71 break; 

72 end 
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74 end 

75 

76 % Compute L2—error and plot the solution 

11 [ 12err , hmax]= dg_error(eoef , mesh , @fexaet ,@fdiff, degree); 

78 

79 % Degree of freedom 

80 dof=size (mesh .Element s,l)*(degree+l)=t:(degree+2)=i<0.5; 

81 

82 fprintf(’ ^DoFs^h_max^L2—error i t \ n ’) 

83 
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84 fprintf (’%7d^%5.3f^%5.3e^%d\n’ 

85 dof , hmax ,12err ,noi); 

86 end % End of function Main_Nonlinear 

87 

88 Define diffusion , advection , and reaction as subfunctions 

89 

90 % Diffusion 

91 function diff = fdiff(x,y) 

92 diff = (10^( —6)). * ones ( size (x )); 

93 end 

94 

95 % Advection 

96 function [advl,adv2] = fadv(x,y) 

97 advl =( 1/sqrt(5))* ones(size(x)); 

98 adv2 =(2/sqrt(5))* ones(size(x)); 

99 end 

100 

101 % Linear reaction 

102 function react = freact(x,y) 

103 react = ones(size(x )); 

104 end 

105 

106 % Nonlinear reaction 

107 function [r,dr] = freact_nonlinear (u) 

108 % Value of the nonlinear reaction term at the current iterate 

109 r = u.^2; 

110 

111 % Value of the derivative of the nonlinear reaction 

112 % term at the current iterate 

113 dr = 2*u; 

114 end 

115 

116 'Wo Define exact solution and force as subfunctions 

117 

118 % Exact solution 

119 function [yex,yex_x,yex_y] = fexact(fdiff ,x,y) 

120 % Evaluate the diffusion function 

121 diff = feval(fdiff ,x,y); 

122 % Exact value 

123 yex =0.5 *(1 — tanh ((2* x—y —0.25)./(sqrt(5*diff )))); 

124 % Eirst derivative wrt x 

125 yex_x = ( — l./(sqrt(5*diff))).*(sech ((2* x—y —0.25)./... 

126 (sqrt(5* diff )))).^2; 

127 % Eirst derivative wrt y 
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128 yex_y = ((0.5)./(sqrt(5*diff))).*(sech((2* x—y —0.25)./... 

129 (sqrt(5*diff )))).^2; 

130 end 

131 

132 % Force function 

133 function source = fsource(fdiff, fadv ,freact ,x,y) 

134 % Evaluate the diffusion function 

135 diff = feval(fdiff ,X,y ); 

136 % Evaluate the advection function 

137 [advl,adv2] = feval(fadv,x, y ); 

138 % Evaluate the linear reaction function 

139 reac = feval(freact ,x,y); 

140 

141 % Exact value 

142 yex=0.5*(l — tanh ((2* x—y —0.25)./(sqrt(5*diff )))); 

143 % First derivative wrt x 

144 yex_x = ( — l./(sqrt(5*diff))).*(sech ((2* x—y —0.25)./... 

145 (sqrt(5* diff )))).^2; 

146 % First derivative wrt y 

147 yex_y = ((0.5)./(sqrt(5*diff))).=t=(sech((2* x—y —0.25)./... 

148 (sqrt(5* diff )))).^2; 

149 % Second derivative wrt x 

150 yex_xx = ((0.8)./diff).* tanh ((2* x—y —0.25)./(sqrt(5*diff ))). * 

151 ( sech ((2* X— y —0.25)./(sqrt(5*diff )))).^2; 

152 % Second derivative wrt y 

153 yex_yy = ((0.2)./ diff).* tanh ((2* x—y — 0.25)./(sqrt(5*diff ))).* 

154 (sech ((2* X— y —0.25)./(sqrt(5*diff )))).^2; 

155 

156 % Force function 

157 source=— diff . * ( yex_xx+yex_yy ) + ( advl . * yex_x + adv2 . * yex_y ) + .. 

158 reac.*yex+yex.^2; 

159 end 

160 
161 

162 Boundary Conditions 

163 

164 % Drichlet Boundary Condition 

165 function DBC=DBCexact ( fdiff , x , y) 

166 % Evaluate the diffusion function 

167 diff = feval(fdiff ,X,y ); 

168 % Drichlet Boundary Condition 

169 DBC=0.5*(1 — tanh ((2* x—y —0.25)./(sqrt(5*diff )))); 

170 end 

171 
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172 % Neumann Boundary Condition 

173 function NC = NBCexact (mesh , fdiff , x , y) 

174 % Neumann Boundary Condition 

175 NC=zeros(size(X )); 

176 end 
111 

178 

179 

180 Set-up parameters function for DG—FEM 

181 

182 function [ penalty , kappa ]= set_parameter (method , degree ) 

183 

184 global Equation ; 

185 

186 Equation . b0 = l; % Superpenalization parameter (In standart b0 = l) 

187 Equation . base =2; % Choose the basis ( 1: monomials, 2:Dubiner Basis) 

188 

189 switeh method 

190 ease 1 

191 %NIPG 

192 Equation . method = 1; 

193 kappa = l; % type of primal method 

194 penalty=l; % penalty parameter 

195 ease 2 

196 

197 Equation . method = 2; 

198 kappa = —1; % type of primal method 

199 penalty =3* degree *( degree + 1); % penalty parameter 

200 ease 3 

201 %UPG 

202 Equation . method = 3; 

203 kappa=0; % type of primal method 

204 penalty =3* degree *( degree + 1); % penalty parameter 

205 

206 end 

207 

208 end 
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